% ML TAU 2013 final project script
% Alt 10:
% Calculate Earth Movers' Distance between each element of X1 and X2,
% then perform SVM on the resultant vector
% Use the Fast EMD algorithm implemented by Ofir Pele, Michael Werman
% Based on ""

%base_path = '/home/itay/TAU/IML/final';

% Add path to the libsvm
%addpath(base_path);
addpath('./libsvm-3.17/matlab');
addpath('./FastEMD');

%
% Initialization
%
close all; clear; clc; tic;
if ~exist('dataforproject.mat','file')
    error('Data file not found in current directory')
end;
% Saving memory. Loading vars on a need to load basis only throughout
load('dataforproject.mat','X1train','X2train','gidtrain','ytrain');
fprintf('Training Data successfully loaded\n');

[ num_features, num_samples ] = size(X1train);
histogram_size = 61;
num_blocks = (num_features / histogram_size);

% EMD parameters
%EMD_threshold= 3;
EMD_extra_mass_penalty = -1;
EMD_cache_file = 'EMD_Mat_SimpleDist.mat';

% Calculate EMD for each two 61 features histograms
% This is a heavy calculation - cache it in file
if ~exist(EMD_cache_file, 'file')

fprintf('Cache file for EMD "%s" not found - performing calculation\n', EMD_cache_file);

EMD_Mat = zeros(num_blocks, num_samples);
for sample = 1:num_samples
	for block = 0:(num_blocks - 1)
		V1 = X1train(histogram_size*block + 1:histogram_size*(block + 1), sample);
		V2 = X2train(histogram_size*block + 1:histogram_size*(block + 1), sample);

		% Create the distance matrix: D_ij = Distance between V1(i) and V2(j)
		%EMD_ground_distance = ones(histogram_size, histogram_size).*EMD_threshold;
		%for i=1:histogram_size
    	%	for j=max([1 i-EMD_threshold+1]):min([histogram_size i+EMD_threshold-1])
        %		EMD_ground_distance(i,j) = abs(V1(i)-V2(j)); 
    	%	end
		%end
		for i=1:histogram_size
			for j=1:histogram_size
				EMD_ground_distance(i, j) = abs(V1(i) - V2(j));
			end
		end

		EMD_Mat(block + 1, sample) = emd_hat_gd_metric_mex(V1, V2, EMD_ground_distance, EMD_extra_mass_penalty);
	end
	fprintf('Done sample %d\n', sample);
end


fprintf('Saving results to cache file "%s"\n', EMD_cache_file);

save(EMD_cache_file, 'EMD_Mat');

else

fprintf('Loading EMD information from "%s"...\n', EMD_cache_file);

load(EMD_cache_file, 'EMD_Mat');

fprintf('EMD data loaded\n');

end

% Smooth the data by applying log
%EMD_Mat = log(EMD_Mat);

c_vals = logspace(0,5,20);

% Optimization - since the data is the same, pre-calcualte
% the test and train data

[ num_emd_features, num_train_entries ] = size(EMD_Mat(:, (gidtrain~=1)));
num_test_entries = size(EMD_Mat(:, (gidtrain==1)), 2);

NormalizedTrainX = zeros(num_train_entries, num_emd_features, 3);
NormalizedTrainX = zeros(num_train_entries, num_emd_features, 3);
y_train = zeros(3, num_train_entries, 1)';
y_test = zeros(3, num_test_entries, 1)';

for k=1:3
	% Training set
	X = EMD_Mat(:, (gidtrain~=k))';
	y_train(:, k) = ytrain(gidtrain~=k);

	% Scale X to [ -1, 1 ]
	Minimums = min(X, [], 1);
	Ranges = max(X, [], 1) - Minimums;
	NormalizedTrainX(:, :, k) = (2 * (X - repmat(Minimums, num_train_entries, 1)) ./ repmat(Ranges, num_train_entries, 1)) - 1; 

	% Testing set
	X = EMD_Mat(:, (gidtrain==k))';	
	y_test(:, k) = ytrain(gidtrain==k);

	% Scale to [ -1, 1 ]
	NormalizedTestX(:, :, k) = (2 * (X - repmat(Minimums, num_test_entries, 1)) ./ repmat(Ranges, num_test_entries, 1)) - 1; 
end

for t=1:3
	for C=1:length(c_vals)
		for d=1:4

Results = zeros(3, 1);
for k=1:3
	
	model = svmtrain(y_train(:, k),NormalizedTrainX(:, :, k), sprintf('-t %d -g 1.0 -c %f -d %d -h 0 -q', t, c_vals(C), d));

	% Predict
	predictedY = svmpredict(y_test(:, k),NormalizedTestX(:, :, k),model);

	Results(k) = sum(y_test(:, k).*predictedY > 0) / size(y_test(:, k),1); 

end

fprintf('Results for kernel type=%d C=%f, d=%d: %f\n', t, C, d, mean(Results));

end
end
end

